APIM - Control - Final Models
Pascal Küng, MSc, Patrick Höhener, MSc, James M. Allen, PhD, Robert Tobias, PhD, Urte Scholz, PhD
2025-03-06
knitr::opts_chunk$set(echo = TRUE)
rm(list = ls())
#cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")
#cmdstanr::set_cmdstan_path(path = "C:/Users/pascku/.cmdstan/cmdstan-2.36.0")
library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)
source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))## [1] 1116# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
do_priorsense = FALSE
get_bayesfactor = TRUE
check_models = TRUE #
if (get_bayesfactor) {
stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'BF', 'Rhat', 'ESS')
} else {
stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS')
}
options(
dplyr.print_max = 100,
brms.backend = 'cmdstan',
brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
brms.file_refit = 'on_change',
#brms.file_refit = 'always',
error = function() {
beepr::beep(sound = 5)
if (shutdown) {
system("shutdown /s /t 180")
quit(save = "no", status = 1)
}
}
, es.use_symbols = TRUE
)
####################### Model parameters #######################
iterations = 12000 # 12'000 per chain to achieve 40'000
warmup = 2000 # 2000
# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'
################################################################
suffix = paste0('_AsPreregisteredAPIM_', as.character(iterations))df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df
df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]Constructing scales Re-coding pusing reshaping data (4field) centering data within and between
Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.1649 0.0000 5.0000 275
Prepare for APIM
Modelling
Self-Reported MVPA
Persuasion
formula <- bf(
pa_sub ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID),
hu = ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_A', dpar = 'hu')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_B', dpar = 'hu')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = hurdle_lognormal()
#)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_sub_persuasion <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::hurdle_lognormal(),
#family = brms::hurdle_negbinomial(),
#family = brms::hurdle_poisson(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 42,
file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_persuasion", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.## Start samplingif (check_models) {
check_brms(pa_sub_persuasion, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_sub_persuasion, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 2.24 [2.14, 2.33] 1.49 0.45
## is_B 1.06 [1.03, 1.10] 1.03 0.95
## is_A:persuasion_self_cw 1.15 [1.12, 1.19] 1.07 0.87
## is_B:persuasion_self_cw 1.14 [1.11, 1.18] 1.07 0.88
## is_A:persuasion_partner_cw 1.04 [1.02, 1.09] 1.02 0.96
## is_B:persuasion_partner_cw 2.71 [2.60, 2.84] 1.65 0.37
## is_A:persuasion_self_cb 2.86 [2.74, 2.99] 1.69 0.35
## Tolerance 95% CI
## [0.43, 0.47]
## [0.91, 0.97]
## [0.84, 0.89]
## [0.85, 0.90]
## [0.92, 0.98]
## [0.35, 0.38]
## [0.33, 0.37]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 94%
## lognormal 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## neg. binomial (zero-infl.) 78%
## beta-binomial 16%
## lognormal 3%
##
## Divergences:## 0 of 40000 iterations ended with a divergence.##
## Tree depth:## 0 of 40000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 7, observations = 3742, p-value =
## 0.84
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.0005344735 0.0029396045
## sample estimates:
## outlier frequency (expected: 0.001574024585783 )
## 0.001870657# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub_persuasion$data$pa_sub[pa_sub_persuasion$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_sub <- summarize_brms(
pa_sub_persuasion,
stats_to_report = stats_to_report,
rope_range = rope_range_continuous,
hu_rope_range = c(-0.18, 0.18),
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Model has no intercept. VIFs may not be sensible.## Sampling priors, please wait...## Warning in summarize_brms(pa_sub_persuasion, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.| exp(Est.)_hu | SE_hu | 95% CI_hu | pd_hu | ROPE_hu | inside ROPE_hu | BF_hu | BF_Evidence_hu | Rhat_hu | Bulk_ESS_hu | Tail_ESS_hu | exp(Est.)_nonzero | SE_nonzero | 95% CI_nonzero | pd_nonzero | ROPE_nonzero | inside ROPE_nonzero | BF_nonzero | BF_Evidence_nonzero | Rhat_nonzero | Bulk_ESS_nonzero | Tail_ESS_nonzero | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.74 | 0.14 | [ 0.51, 1.07] | 0.944 | [0.84, 1.20] | 0.257 | 0.260 | Moderate Evidence for Null | 1.000 | 17468 | 24772 | 46.58*** | 3.50 | [40.07, 53.85] | 1.000 | [0.92, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 13642 | 22685 |
| is_A:day | 0.97 | 0.18 | [ 0.67, 1.40] | 0.570 | [0.84, 1.20] | 0.660 | 0.095 | Strong Evidence for Null | 1.000 | 44612 | 29960 | 1.09 | 0.09 | [ 0.92, 1.29] | 0.845 | [0.92, 1.08] | 0.436 | 0.012 | Very Strong Evidence for Null | 1.000 | 46827 | 29347 |
| is_A:persuasion_partner_cb | 0.75 | 0.30 | [ 0.34, 1.64] | 0.767 | [0.84, 1.20] | 0.271 | 0.254 | Moderate Evidence for Null | 1.000 | 14088 | 21142 | 1.02 | 0.17 | [ 0.74, 1.42] | 0.561 | [0.92, 1.08] | 0.365 | 0.018 | Very Strong Evidence for Null | 1.000 | 11377 | 19396 |
| is_A:persuasion_partner_cw | 1.58*** | 0.16 | [ 1.30, 1.98] | 1.000 | [0.84, 1.20] | 0.004 | >100 | Overwhelming Evidence | 1.000 | 23821 | 26789 | 1.02 | 0.03 | [ 0.96, 1.08] | 0.712 | [0.92, 1.08] | 0.980 | 0.008 | Very Strong Evidence for Null | 1.000 | 20463 | 25811 |
| is_A:persuasion_self_cb | 0.81 | 0.36 | [ 0.33, 1.97] | 0.683 | [0.84, 1.20] | 0.281 | 0.251 | Moderate Evidence for Null | 1.000 | 13340 | 20895 | 1.21 | 0.21 | [ 0.84, 1.72] | 0.853 | [0.92, 1.08] | 0.199 | 0.029 | Very Strong Evidence for Null | 1.001 | 10250 | 19030 |
| is_A:persuasion_self_cw | 1.67*** | 0.16 | [ 1.41, 2.10] | 1.000 | [0.84, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 22586 | 23823 | 1.01 | 0.04 | [ 0.94, 1.10] | 0.643 | [0.92, 1.08] | 0.936 | 0.011 | Very Strong Evidence for Null | 1.000 | 15189 | 22445 |
| is_B | 0.92 | 0.18 | [ 0.63, 1.37] | 0.656 | [0.84, 1.20] | 0.597 | 0.084 | Strong Evidence for Null | 1.000 | 16516 | 22811 | 49.76*** | 3.97 | [42.49, 58.38] | 1.000 | [0.92, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 11733 | 19760 |
| is_B:day | 0.87 | 0.16 | [ 0.61, 1.25] | 0.777 | [0.84, 1.20] | 0.543 | 0.123 | Moderate Evidence for Null | 1.000 | 41627 | 29767 | 0.88 | 0.07 | [ 0.75, 1.03] | 0.943 | [0.92, 1.08] | 0.257 | 0.024 | Very Strong Evidence for Null | 1.000 | 50533 | 30069 |
| is_B:persuasion_partner_cb | 1.01 | 0.48 | [ 0.39, 2.62] | 0.510 | [0.84, 1.20] | 0.293 | 0.235 | Moderate Evidence for Null | 1.000 | 12851 | 19993 | 1.03 | 0.19 | [ 0.72, 1.48] | 0.568 | [0.92, 1.08] | 0.330 | 0.017 | Very Strong Evidence for Null | 1.000 | 11737 | 19703 |
| is_B:persuasion_partner_cw | 1.49*** | 0.14 | [ 1.24, 1.85] | 1.000 | [0.84, 1.20] | 0.010 | 86.845 | Very Strong Evidence | 1.000 | 26845 | 24277 | 1.03 | 0.03 | [ 0.97, 1.10] | 0.858 | [0.92, 1.08] | 0.935 | 0.014 | Very Strong Evidence for Null | 1.000 | 20439 | 27855 |
| is_B:persuasion_self_cb | 0.70 | 0.29 | [ 0.30, 1.62] | 0.805 | [0.84, 1.20] | 0.233 | 0.304 | Weak Evidence for Null | 1.000 | 13544 | 21173 | 1.13 | 0.18 | [ 0.82, 1.57] | 0.767 | [0.92, 1.08] | 0.289 | 0.023 | Very Strong Evidence for Null | 1.000 | 12920 | 20807 |
| is_B:persuasion_self_cw | 1.68*** | 0.17 | [ 1.39, 2.10] | 1.000 | [0.84, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 22677 | 24944 | 1.04 | 0.03 | [ 0.98, 1.10] | 0.889 | [0.92, 1.08] | 0.930 | 0.016 | Very Strong Evidence for Null | 1.000 | 17302 | 24439 |
| sd(is_A) | 0.94 | 0.12 | [0.73, 1.23] | NA | NA | NA | NA | NA | 1.000 | 20183 | 25071 | 0.33 | 0.05 | [0.25, 0.44] | NA | NA | NA | NA | NA | 1.000 | 13315 | 23621 |
| sd(is_A:persuasion_partner_cw) | 0.41 | 0.12 | [0.21, 0.69] | NA | NA | NA | NA | NA | 1.000 | 18656 | 21846 | 0.06 | 0.03 | [0.01, 0.13] | NA | NA | NA | NA | NA | 1.000 | 12477 | 10932 |
| sd(is_A:persuasion_self_cw) | 0.29 | 0.13 | [0.05, 0.58] | NA | NA | NA | NA | NA | 1.000 | 9601 | 8514 | 0.17 | 0.03 | [0.11, 0.24] | NA | NA | NA | NA | NA | 1.000 | 20343 | 28866 |
| sd(is_B) | 1.01 | 0.13 | [0.79, 1.32] | NA | NA | NA | NA | NA | 1.000 | 18297 | 24221 | 0.39 | 0.05 | [0.30, 0.50] | NA | NA | NA | NA | NA | 1.000 | 13771 | 21334 |
| sd(is_B:persuasion_partner_cw) | 0.33 | 0.14 | [0.06, 0.63] | NA | NA | NA | NA | NA | 1.000 | 10537 | 8117 | 0.10 | 0.03 | [0.05, 0.16] | NA | NA | NA | NA | NA | 1.000 | 20368 | 18826 |
| sd(is_B:persuasion_self_cw) | 0.38 | 0.10 | [0.19, 0.62] | NA | NA | NA | NA | NA | 1.000 | 21344 | 22705 | 0.09 | 0.03 | [0.05, 0.15] | NA | NA | NA | NA | NA | 1.000 | 18840 | 17891 |
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 0.67 | 0.01 | [0.65, 0.70] | NA | NA | NA | NA | NA | 1.000 | 47834 | 30606 |
Hurdle part of the model on the left, non-zero part towards the right side of the table
Pressure
formula <- bf(
pa_sub ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID),
hu = ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_A', dpar = 'hu')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_B', dpar = 'hu')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = hurdle_lognormal()
#)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_sub_pressure <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::hurdle_lognormal(),
#family = brms::hurdle_negbinomial(),
#family = brms::hurdle_poisson(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 42,
file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_pressure", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.## Start samplingif (check_models) {
check_brms(pa_sub_pressure, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_sub_pressure, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 2.23 [2.14, 2.33] 1.49 0.45
## is_B 1.05 [1.02, 1.09] 1.02 0.96
## is_A:pressure_self_cw 1.01 [1.00, 1.14] 1.01 0.99
## is_B:pressure_self_cw 1.01 [1.00, 1.12] 1.01 0.99
## is_A:pressure_partner_cw 1.01 [1.00, 1.15] 1.01 0.99
## is_B:pressure_partner_cw 3.49 [3.33, 3.66] 1.87 0.29
## is_A:pressure_self_cb 4.30 [4.10, 4.52] 2.07 0.23
## Tolerance 95% CI
## [0.43, 0.47]
## [0.92, 0.98]
## [0.88, 1.00]
## [0.89, 1.00]
## [0.87, 1.00]
## [0.27, 0.30]
## [0.22, 0.24]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 94%
## lognormal 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## neg. binomial (zero-infl.) 78%
## beta-binomial 16%
## lognormal 3%
##
## Divergences:## 0 of 40000 iterations ended with a divergence.##
## Tree depth:## 0 of 40000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 11 observations with a pareto_k > 0.7 in model 'model'.
## We recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 6, observations = 3736, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.0005353319 0.0030848501
## sample estimates:
## outlier frequency (expected: 0.00157922912205567 )
## 0.001605996# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub_pressure$data$pa_sub[pa_sub_pressure$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_sub <- summarize_brms(
pa_sub_pressure,
stats_to_report = stats_to_report,
rope_range = rope_range_continuous,
hu_rope_range = c(-0.18, 0.18),
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Model has no intercept. VIFs may not be sensible.## Sampling priors, please wait...## Warning in summarize_brms(pa_sub_pressure, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.| exp(Est.)_hu | SE_hu | 95% CI_hu | pd_hu | ROPE_hu | inside ROPE_hu | BF_hu | BF_Evidence_hu | Rhat_hu | Bulk_ESS_hu | Tail_ESS_hu | exp(Est.)_nonzero | SE_nonzero | 95% CI_nonzero | pd_nonzero | ROPE_nonzero | inside ROPE_nonzero | BF_nonzero | BF_Evidence_nonzero | Rhat_nonzero | Bulk_ESS_nonzero | Tail_ESS_nonzero | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.87 | 0.15 | [ 0.61, 1.23] | 0.783 | [0.84, 1.20] | 0.555 | 0.096 | Strong Evidence for Null | 1.000 | 14549 | 20069 | 48.25*** | 3.51 | [41.79, 55.67] | 1.000 | [0.92, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 10797 | 19767 |
| is_A:day | 0.67* | 0.12 | [ 0.48, 0.95] | 0.989 | [0.84, 1.20] | 0.107 | 1.172 | Weak Evidence | 1.000 | 37555 | 29688 | 1.09 | 0.09 | [ 0.92, 1.28] | 0.835 | [0.92, 1.08] | 0.452 | 0.011 | Very Strong Evidence for Null | 1.000 | 48844 | 30086 |
| is_A:pressure_partner_cb | 0.63 | 0.39 | [ 0.18, 2.10] | 0.772 | [0.84, 1.20] | 0.174 | 0.412 | Weak Evidence for Null | 1.000 | 14001 | 22562 | 0.94 | 0.26 | [ 0.54, 1.64] | 0.582 | [0.92, 1.08] | 0.218 | 0.016 | Very Strong Evidence for Null | 1.000 | 12801 | 20013 |
| is_A:pressure_partner_cw | 2.44*** | 0.57 | [ 1.57, 4.32] | 1.000 | [0.84, 1.20] | 0.001 | >100 | Overwhelming Evidence | 1.000 | 27317 | 18821 | 0.93 | 0.05 | [ 0.84, 1.05] | 0.880 | [0.92, 1.08] | 0.567 | 0.012 | Very Strong Evidence for Null | 1.000 | 29892 | 25384 |
| is_A:pressure_self_cb | 0.57 | 0.44 | [ 0.13, 2.62] | 0.766 | [0.84, 1.20] | 0.143 | 0.502 | Weak Evidence for Null | 1.000 | 14022 | 21390 | 1.56 | 0.53 | [ 0.80, 3.08] | 0.908 | [0.92, 1.08] | 0.075 | 0.038 | Strong Evidence for Null | 1.000 | 14697 | 21985 |
| is_A:pressure_self_cw | 1.86* | 0.62 | [ 1.02, 4.68] | 0.978 | [0.84, 1.20] | 0.060 | 1.518 | Weak Evidence | 1.000 | 19579 | 18590 | 0.85* | 0.06 | [ 0.72, 1.00] | 0.978 | [0.92, 1.08] | 0.136 | 0.082 | Strong Evidence for Null | 1.000 | 23058 | 21388 |
| is_B | 1.14 | 0.21 | [ 0.80, 1.63] | 0.763 | [0.84, 1.20] | 0.560 | 0.090 | Strong Evidence for Null | 1.000 | 13997 | 21736 | 52.20*** | 4.30 | [44.46, 61.41] | 1.000 | [0.92, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 9866 | 16264 |
| is_B:day | 0.56*** | 0.10 | [ 0.40, 0.79] | 1.000 | [0.84, 1.20] | 0.011 | 20.363 | Strong Evidence | 1.000 | 36258 | 29993 | 0.83* | 0.07 | [ 0.71, 0.97] | 0.991 | [0.92, 1.08] | 0.085 | 0.109 | Moderate Evidence for Null | 1.000 | 48596 | 30206 |
| is_B:pressure_partner_cb | 1.13 | 0.88 | [ 0.24, 5.28] | 0.562 | [0.84, 1.20] | 0.183 | 0.393 | Weak Evidence for Null | 1.000 | 12957 | 20347 | 1.10 | 0.43 | [ 0.49, 2.47] | 0.593 | [0.92, 1.08] | 0.155 | 0.018 | Very Strong Evidence for Null | 1.001 | 12994 | 19397 |
| is_B:pressure_partner_cw | 3.28* | 2.01 | [ 1.24, 16.35] | 0.990 | [0.84, 1.20] | 0.017 | 4.367 | Moderate Evidence | 1.000 | 17196 | 17609 | 0.99 | 0.07 | [ 0.82, 1.13] | 0.575 | [0.92, 1.08] | 0.709 | 0.008 | Very Strong Evidence for Null | 1.000 | 21016 | 19395 |
| is_B:pressure_self_cb | 0.48 | 0.30 | [ 0.14, 1.67] | 0.882 | [0.84, 1.20] | 0.111 | 0.645 | Weak Evidence for Null | 1.000 | 12815 | 19815 | 1.08 | 0.34 | [ 0.57, 2.05] | 0.593 | [0.92, 1.08] | 0.194 | 0.019 | Very Strong Evidence for Null | 1.001 | 12751 | 19790 |
| is_B:pressure_self_cw | 1.29 | 0.34 | [ 0.71, 2.66] | 0.832 | [0.84, 1.20] | 0.326 | 0.227 | Moderate Evidence for Null | 1.000 | 20759 | 17702 | 0.93 | 0.06 | [ 0.82, 1.08] | 0.842 | [0.92, 1.08] | 0.545 | 0.012 | Very Strong Evidence for Null | 1.000 | 23751 | 24390 |
| sd(is_A) | 0.88 | 0.12 | [0.68, 1.15] | NA | NA | NA | NA | NA | 1.000 | 17349 | 24619 | 0.32 | 0.05 | [0.24, 0.44] | NA | NA | NA | NA | NA | 1.000 | 12153 | 20889 |
| sd(is_A:pressure_partner_cw) | 0.31 | 0.28 | [0.01, 1.20] | NA | NA | NA | NA | NA | 1.000 | 15331 | 16422 | 0.06 | 0.05 | [0.00, 0.21] | NA | NA | NA | NA | NA | 1.000 | 17943 | 17235 |
| sd(is_A:pressure_self_cw) | 0.72 | 0.53 | [0.06, 2.40] | NA | NA | NA | NA | NA | 1.000 | 11421 | 12927 | 0.10 | 0.09 | [0.00, 0.38] | NA | NA | NA | NA | NA | 1.000 | 12722 | 15622 |
| sd(is_B) | 0.91 | 0.12 | [0.70, 1.19] | NA | NA | NA | NA | NA | 1.000 | 18280 | 24673 | 0.40 | 0.05 | [0.31, 0.53] | NA | NA | NA | NA | NA | 1.000 | 11041 | 18555 |
| sd(is_B:pressure_partner_cw) | 1.65 | 0.85 | [0.31, 3.83] | NA | NA | NA | NA | NA | 1.000 | 12805 | 11132 | 0.09 | 0.08 | [0.00, 0.37] | NA | NA | NA | NA | NA | 1.000 | 13239 | 14438 |
| sd(is_B:pressure_self_cw) | 0.74 | 0.49 | [0.06, 2.12] | NA | NA | NA | NA | NA | 1.000 | 10281 | 10816 | 0.09 | 0.07 | [0.00, 0.28] | NA | NA | NA | NA | NA | 1.000 | 14609 | 13572 |
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 0.69 | 0.01 | [0.67, 0.72] | NA | NA | NA | NA | NA | 1.000 | 43470 | 29205 |
Hurdle part of the model on the left, non-zero part towards the right side of the table
Pushing
formula <- bf(
pa_sub ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
barriers_self_cw + barriers_self_cb +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID),
hu = ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_A', dpar = 'hu')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_B', dpar = 'hu')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = hurdle_lognormal()
#)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_sub_pushing <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::hurdle_lognormal(),
#family = brms::hurdle_negbinomial(),
#family = brms::hurdle_poisson(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 42,
file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_pushing", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.## Start samplingif (check_models) {
check_brms(pa_sub_pushing, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_sub_pushing, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.75 [1.68, 1.84] 1.32 0.57
## is_B 1.08 [1.05, 1.12] 1.04 0.93
## barriers_self_cw 1.04 [1.02, 1.10] 1.02 0.96
## barriers_self_cb 1.04 [1.01, 1.09] 1.02 0.97
## is_A:pushing_self_cw 1.07 [1.04, 1.12] 1.04 0.93
## is_B:pushing_self_cw 1.22 [1.18, 1.28] 1.11 0.82
## is_A:pushing_partner_cw 1.24 [1.20, 1.30] 1.12 0.80
## is_B:pushing_partner_cw 1.19 [1.15, 1.24] 1.09 0.84
## is_A:pushing_self_cb 1.27 [1.23, 1.33] 1.13 0.79
## Tolerance 95% CI
## [0.54, 0.60]
## [0.89, 0.95]
## [0.91, 0.98]
## [0.92, 0.99]
## [0.89, 0.96]
## [0.78, 0.85]
## [0.77, 0.83]
## [0.81, 0.87]
## [0.75, 0.82]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 91%
## lognormal 9%
##
## Predicted Distribution of Response
##
## Distribution Probability
## neg. binomial (zero-infl.) 84%
## beta-binomial 9%
## lognormal 6%
##
## Divergences:## 0 of 40000 iterations ended with a divergence.##
## Tree depth:## 0 of 40000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 12 observations with a pareto_k > 0.7 in model 'model'.
## We recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 5, observations = 1776, p-value =
## 0.36
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.00423705
## sample estimates:
## outlier frequency (expected: 0.00161036036036036 )
## 0.002815315# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub_pushing$data$pa_sub[pa_sub_pushing$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_sub <- summarize_brms(
pa_sub_pushing,
stats_to_report = stats_to_report,
rope_range = rope_range_continuous,
hu_rope_range = c(-0.18, 0.18),
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Model has no intercept. VIFs may not be sensible.## Sampling priors, please wait...## Warning in summarize_brms(pa_sub_pushing, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.| exp(Est.)_hu | SE_hu | 95% CI_hu | pd_hu | ROPE_hu | inside ROPE_hu | BF_hu | BF_Evidence_hu | Rhat_hu | Bulk_ESS_hu | Tail_ESS_hu | exp(Est.)_nonzero | SE_nonzero | 95% CI_nonzero | pd_nonzero | ROPE_nonzero | inside ROPE_nonzero | BF_nonzero | BF_Evidence_nonzero | Rhat_nonzero | Bulk_ESS_nonzero | Tail_ESS_nonzero | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| barriers_self_cb | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 0.89 | 0.13 | [ 0.68, 1.18] | 0.791 | [0.93, 1.08] | 0.305 | 0.026 | Very Strong Evidence for Null | 1.000 | 15389 | 22099 |
| barriers_self_cw | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 1.06 | 0.04 | [ 0.99, 1.13] | 0.943 | [0.93, 1.08] | 0.739 | 0.032 | Strong Evidence for Null | 1.000 | 68293 | 29503 |
| is_A | 2.63*** | 0.70 | [ 1.56, 4.51] | 1.000 | [0.84, 1.20] | 0.002 | 48.662 | Very Strong Evidence | 1.000 | 29400 | 28949 | 53.39*** | 4.54 | [45.24, 63.19] | 1.000 | [0.93, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 18788 | 26280 |
| is_A:day | 0.82 | 0.24 | [ 0.46, 1.45] | 0.751 | [0.84, 1.20] | 0.369 | 0.187 | Moderate Evidence for Null | 1.000 | 58113 | 28964 | 1.08 | 0.10 | [ 0.90, 1.30] | 0.793 | [0.93, 1.08] | 0.441 | 0.011 | Very Strong Evidence for Null | 1.000 | 66370 | 30302 |
| is_A:pushing_partner_cb | 0.91 | 0.88 | [ 0.13, 6.16] | 0.540 | [0.84, 1.20] | 0.147 | 0.491 | Weak Evidence for Null | 1.000 | 30887 | 30885 | 1.10 | 0.40 | [ 0.52, 2.26] | 0.597 | [0.93, 1.08] | 0.159 | 0.020 | Very Strong Evidence for Null | 1.000 | 15843 | 23109 |
| is_A:pushing_partner_cw | 1.95** | 0.45 | [ 1.30, 3.48] | 0.999 | [0.84, 1.20] | 0.009 | 14.509 | Strong Evidence | 1.000 | 32690 | 26682 | 0.96 | 0.03 | [ 0.89, 1.03] | 0.886 | [0.93, 1.08] | 0.837 | 0.017 | Very Strong Evidence for Null | 1.000 | 47103 | 27861 |
| is_A:pushing_self_cb | 0.39 | 0.35 | [ 0.07, 2.31] | 0.852 | [0.84, 1.20] | 0.093 | 0.769 | Weak Evidence for Null | 1.000 | 28601 | 30135 | 1.29 | 0.42 | [ 0.67, 2.49] | 0.785 | [0.93, 1.08] | 0.136 | 0.025 | Very Strong Evidence for Null | 1.000 | 18803 | 26620 |
| is_A:pushing_self_cw | 1.91 | 0.76 | [ 0.96, 5.60] | 0.967 | [0.84, 1.20] | 0.081 | 0.998 | Weak Evidence for Null | 1.000 | 25693 | 21815 | 1.00 | 0.06 | [ 0.89, 1.12] | 0.523 | [0.93, 1.08] | 0.814 | 0.012 | Very Strong Evidence for Null | 1.000 | 24382 | 27449 |
| is_B | 4.02*** | 1.18 | [ 2.24, 7.28] | 1.000 | [0.84, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 26916 | 28674 | 52.54*** | 4.45 | [44.44, 62.18] | 1.000 | [0.93, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 20116 | 25664 |
| is_B:day | 0.65 | 0.19 | [ 0.36, 1.15] | 0.929 | [0.84, 1.20] | 0.177 | 0.427 | Weak Evidence for Null | 1.000 | 63463 | 30099 | 0.94 | 0.09 | [ 0.79, 1.12] | 0.746 | [0.93, 1.08] | 0.507 | 0.010 | Very Strong Evidence for Null | 1.000 | 74905 | 30566 |
| is_B:pushing_partner_cb | 0.81 | 0.81 | [ 0.11, 5.67] | 0.582 | [0.84, 1.20] | 0.140 | 0.508 | Weak Evidence for Null | 1.000 | 29944 | 29697 | 1.04 | 0.34 | [ 0.53, 2.01] | 0.543 | [0.93, 1.08] | 0.187 | 0.020 | Very Strong Evidence for Null | 1.000 | 19801 | 25578 |
| is_B:pushing_partner_cw | 2.28* | 1.11 | [ 1.04, 8.55] | 0.979 | [0.84, 1.20] | 0.044 | 1.868 | Weak Evidence | 1.000 | 21839 | 21728 | 1.02 | 0.06 | [ 0.91, 1.15] | 0.658 | [0.93, 1.08] | 0.797 | 0.012 | Very Strong Evidence for Null | 1.000 | 27666 | 30385 |
| is_B:pushing_self_cb | 0.41 | 0.42 | [ 0.05, 3.18] | 0.810 | [0.84, 1.20] | 0.096 | 0.757 | Weak Evidence for Null | 1.000 | 31398 | 29450 | 1.50 | 0.55 | [ 0.71, 3.13] | 0.859 | [0.93, 1.08] | 0.092 | 0.038 | Strong Evidence for Null | 1.000 | 17644 | 25436 |
| is_B:pushing_self_cw | 1.16 | 0.25 | [ 0.73, 1.91] | 0.753 | [0.84, 1.20] | 0.491 | 0.137 | Moderate Evidence for Null | 1.000 | 29855 | 26726 | 0.95 | 0.03 | [ 0.88, 1.01] | 0.939 | [0.93, 1.08] | 0.726 | 0.030 | Very Strong Evidence for Null | 1.000 | 43301 | 31543 |
| sd(is_A) | 1.18 | 0.20 | [0.85, 1.64] | NA | NA | NA | NA | NA | 1.000 | 27290 | 31419 | 0.36 | 0.06 | [0.26, 0.50] | NA | NA | NA | NA | NA | 1.000 | 16350 | 24786 |
| sd(is_A:pushing_partner_cw) | 0.65 | 0.33 | [0.09, 1.45] | NA | NA | NA | NA | NA | 1.000 | 14902 | 13091 | 0.04 | 0.04 | [0.00, 0.14] | NA | NA | NA | NA | NA | 1.000 | 20863 | 21692 |
| sd(is_A:pushing_self_cw) | 1.51 | 0.59 | [0.49, 3.01] | NA | NA | NA | NA | NA | 1.000 | 13003 | 10459 | 0.18 | 0.06 | [0.05, 0.32] | NA | NA | NA | NA | NA | 1.000 | 12931 | 9839 |
| sd(is_B) | 1.37 | 0.23 | [0.99, 1.92] | NA | NA | NA | NA | NA | 1.000 | 25542 | 30110 | 0.38 | 0.06 | [0.28, 0.52] | NA | NA | NA | NA | NA | 1.000 | 16644 | 25124 |
| sd(is_B:pushing_partner_cw) | 1.72 | 0.82 | [0.24, 3.65] | NA | NA | NA | NA | NA | 1.000 | 10498 | 8114 | 0.17 | 0.07 | [0.03, 0.32] | NA | NA | NA | NA | NA | 1.000 | 12628 | 10344 |
| sd(is_B:pushing_self_cw) | 0.85 | 0.33 | [0.29, 1.72] | NA | NA | NA | NA | NA | 1.000 | 14504 | 13114 | 0.05 | 0.04 | [0.00, 0.15] | NA | NA | NA | NA | NA | 1.000 | 19509 | 21091 |
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 0.67 | 0.01 | [0.64, 0.70] | NA | NA | NA | NA | NA | 1.000 | 57411 | 29137 |
Hurdle part of the model on the left, non-zero part towards the right side of the table
Device Based MVPA
## [1] 5.75 971.25
Persuasion
formula <- bf(
pa_obj ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
is_A:weartime_self_cw + is_B:weartime_self_cw +
is_A:weartime_self_cb + is_B:weartime_self_cb +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = lognormal()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_obj_log_persuasion <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = lognormal(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_persuasion", suffix))
)## Warning: Rows containing NAs were excluded from the model.if (check_models) {
check_brms(pa_obj_log_persuasion, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_obj_log_persuasion, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.84 [1.76, 1.93] 1.36 0.54
## is_B 1.08 [1.05, 1.12] 1.04 0.93
## is_A:persuasion_self_cw 1.13 [1.10, 1.18] 1.06 0.88
## is_B:persuasion_self_cw 1.09 [1.06, 1.13] 1.04 0.92
## is_A:persuasion_partner_cw 1.04 [1.02, 1.09] 1.02 0.96
## is_B:persuasion_partner_cw 3.53 [3.35, 3.73] 1.88 0.28
## is_A:persuasion_self_cb 2.94 [2.79, 3.10] 1.71 0.34
## is_B:persuasion_self_cb 3.69 [3.50, 3.90] 1.92 0.27
## is_A:persuasion_partner_cb 2.95 [2.80, 3.11] 1.72 0.34
## Tolerance 95% CI
## [0.52, 0.57]
## [0.89, 0.95]
## [0.85, 0.91]
## [0.88, 0.95]
## [0.91, 0.98]
## [0.27, 0.30]
## [0.32, 0.36]
## [0.26, 0.29]
## [0.32, 0.36]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 88%
## lognormal 9%
## weibull 3%
##
## Predicted Distribution of Response
##
## Distribution Probability
## lognormal 72%
## tweedie 16%
## neg. binomial (zero-infl.) 9%
##
## Divergences:## 0 of 40000 iterations ended with a divergence.##
## Tree depth:## 0 of 40000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 33, observations = 3343, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.001338618 0.005085253
## sample estimates:
## outlier frequency (expected: 0.00303021238408615 )
## 0.009871373# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log_persuasion$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_obj <- summarize_brms(
pa_obj_log_persuasion,
stats_to_report = stats_to_report,
rope_range = rope_range_log,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.## Sampling priors, please wait...## Warning in summarize_brms(pa_obj_log_persuasion, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 101.53*** | 6.28 | [ 89.80, 114.88] | 1.000 | [0.94, 1.07] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 9441 | 16504 |
| is_B | 135.95*** | 7.32 | [122.13, 151.50] | 1.000 | [0.94, 1.07] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 13358 | 20795 |
| is_A:persuasion_self_cw | 1.05* | 0.02 | [ 1.01, 1.09] | 0.989 | [0.94, 1.07] | 0.839 | 0.083 | Strong Evidence for Null | 1.000 | 36995 | 32045 |
| is_B:persuasion_self_cw | 1.01 | 0.02 | [ 0.97, 1.05] | 0.659 | [0.94, 1.07] | 0.998 | 0.005 | Very Strong Evidence for Null | 1.000 | 35355 | 33024 |
| is_A:persuasion_partner_cw | 1.02 | 0.02 | [ 0.98, 1.06] | 0.792 | [0.94, 1.07] | 0.988 | 0.007 | Very Strong Evidence for Null | 1.000 | 36380 | 30708 |
| is_B:persuasion_partner_cw | 1.03 | 0.02 | [ 0.99, 1.07] | 0.909 | [0.94, 1.07] | 0.969 | 0.014 | Very Strong Evidence for Null | 1.000 | 40799 | 30254 |
| is_A:persuasion_self_cb | 1.25 | 0.21 | [ 0.89, 1.76] | 0.901 | [0.94, 1.07] | 0.128 | 0.033 | Strong Evidence for Null | 1.000 | 10402 | 17765 |
| is_B:persuasion_self_cb | 0.89 | 0.11 | [ 0.70, 1.14] | 0.824 | [0.94, 1.07] | 0.263 | 0.019 | Very Strong Evidence for Null | 1.000 | 14430 | 21279 |
| is_A:persuasion_partner_cb | 0.89 | 0.13 | [ 0.66, 1.20] | 0.787 | [0.94, 1.07] | 0.249 | 0.019 | Very Strong Evidence for Null | 1.000 | 11805 | 19041 |
| is_B:persuasion_partner_cb | 1.29 | 0.19 | [ 0.96, 1.73] | 0.957 | [0.94, 1.07] | 0.081 | 0.062 | Strong Evidence for Null | 1.000 | 13364 | 21848 |
| is_A:day | 0.98 | 0.05 | [ 0.90, 1.08] | 0.628 | [0.94, 1.07] | 0.798 | 0.004 | Very Strong Evidence for Null | 1.000 | 87762 | 28942 |
| is_B:day | 0.95 | 0.04 | [ 0.87, 1.04] | 0.859 | [0.94, 1.07] | 0.612 | 0.007 | Very Strong Evidence for Null | 1.000 | 93845 | 28595 |
| is_A:weartime_self_cw | 1.00*** | 0.00 | [ 1.00, 1.00] | 1.000 | [0.94, 1.07] | 1.000 | 11.432 | Strong Evidence | 1.000 | 79311 | 30606 |
| is_B:weartime_self_cw | 1.00* | 0.00 | [ 1.00, 1.00] | 0.983 | [0.94, 1.07] | 1.000 | 0.039 | Strong Evidence for Null | 1.000 | 80060 | 29963 |
| is_A:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.712 | [0.94, 1.07] | 1.000 | 0.014 | Very Strong Evidence for Null | 1.000 | 12377 | 18753 |
| is_B:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.584 | [0.94, 1.07] | 1.000 | 0.011 | Very Strong Evidence for Null | 1.000 | 15565 | 23274 |
| sd(is_A) | 0.34 | 0.04 | [0.26, 0.44] | NA | NA | NA | NA | NA | 1.000 | 12417 | 22026 |
| sd(is_B) | 0.28 | 0.04 | [0.22, 0.37] | NA | NA | NA | NA | NA | 1.000 | 14817 | 22063 |
| sd(is_A:persuasion_self_cw) | 0.06 | 0.02 | [0.02, 0.11] | NA | NA | NA | NA | NA | 1.000 | 14879 | 11539 |
| sd(is_B:persuasion_self_cw) | 0.06 | 0.02 | [0.02, 0.11] | NA | NA | NA | NA | NA | 1.000 | 17283 | 14227 |
| sd(is_A:persuasion_partner_cw) | 0.07 | 0.02 | [0.03, 0.12] | NA | NA | NA | NA | NA | 1.000 | 15202 | 15518 |
| sd(is_B:persuasion_partner_cw) | 0.07 | 0.03 | [0.01, 0.13] | NA | NA | NA | NA | NA | 1.000 | 11803 | 13070 |
| sigma | 0.55 | 0.01 | [0.53, 0.56] | NA | NA | NA | NA | NA | 1.000 | 73910 | 28826 |
Pressure
formula <- bf(
pa_obj ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
is_A:weartime_self_cw + is_B:weartime_self_cw +
is_A:weartime_self_cb + is_B:weartime_self_cb +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = lognormal()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_obj_log_pressure <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = lognormal(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_pressure", suffix))
)## Warning: Rows containing NAs were excluded from the model.if (check_models) {
check_brms(pa_obj_log_pressure, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_obj_log_pressure, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model
## without interaction terms.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.86 [ 1.78, 1.95] 1.36 0.54
## is_B 1.02 [ 1.01, 1.09] 1.01 0.98
## is_A:pressure_self_cw 1.01 [ 1.00, 1.13] 1.01 0.99
## is_B:pressure_self_cw 1.02 [ 1.00, 1.11] 1.01 0.98
## is_A:pressure_partner_cw 1.04 [ 1.01, 1.09] 1.02 0.96
## Tolerance 95% CI
## [0.51, 0.56]
## [0.91, 0.99]
## [0.88, 1.00]
## [0.90, 1.00]
## [0.92, 0.99]
##
## Moderate Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A:pressure_self_cb 8.49 [ 8.01, 9.01] 2.91 0.12
## is_A:pressure_partner_cb 8.81 [ 8.30, 9.34] 2.97 0.11
## Tolerance 95% CI
## [0.11, 0.12]
## [0.11, 0.12]
##
## High Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_B:pressure_partner_cw 10.71 [10.10, 11.37] 3.27 0.09
## is_B:pressure_self_cb 10.70 [10.08, 11.35] 3.27 0.09
## Tolerance 95% CI
## [0.09, 0.10]
## [0.09, 0.10]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 88%
## lognormal 9%
## weibull 3%
##
## Predicted Distribution of Response
##
## Distribution Probability
## lognormal 72%
## tweedie 16%
## neg. binomial (zero-infl.) 9%
##
## Divergences:## 2 of 40000 iterations ended with a divergence (0.005%).
## Try increasing 'adapt_delta' to remove the divergences.##
## Tree depth:## 0 of 40000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 34, observations = 3337, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.001640695 0.004794726
## sample estimates:
## outlier frequency (expected: 0.00308660473479173 )
## 0.01018879# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log_pressure$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_obj <- summarize_brms(
pa_obj_log_pressure,
stats_to_report = stats_to_report,
rope_range = rope_range_log,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Warning: There were 2 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Model has no intercept. VIFs may not be sensible.## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model
## without interaction terms.## Warning in compute_rope(range = rope_range): Collinearity detected. Some
## VIFs are > 10. This may invalidate ROPE inferences!Check for Multicollinearity
Low Correlation
Term VIF VIF 95% CI Increased SE Tolerance
is_A 1.86 [ 1.78, 1.95] 1.36 0.54
is_B 1.02 [ 1.01, 1.09] 1.01 0.98
is_A:pressure_self_cw 1.01 [ 1.00, 1.13] 1.01 0.99
is_B:pressure_self_cw 1.02 [ 1.00, 1.11] 1.01 0.98is_A:pressure_partner_cw 1.04 [ 1.01, 1.09] 1.02 0.96 Tolerance 95% CI [0.51, 0.56] [0.91, 0.99] [0.88, 1.00] [0.90, 1.00] [0.92, 0.99]
Moderate Correlation
Term VIF VIF 95% CI Increased SE Tolerance
is_A:pressure_self_cb 8.49 [ 8.01, 9.01] 2.91 0.12is_A:pressure_partner_cb 8.81 [ 8.30, 9.34] 2.97 0.11 Tolerance 95% CI [0.11, 0.12] [0.11, 0.12]
High Correlation
Term VIF VIF 95% CI Increased SE Toleranceis_B:pressure_partner_cw 10.71 [10.10, 11.37] 3.27 0.09 is_B:pressure_self_cb 10.70 [10.08, 11.35] 3.27 0.09 Tolerance 95% CI [0.09, 0.10] [0.09, 0.10]
## Sampling priors, please wait...## Warning: Bayes factors might not be precise.
## For precise Bayes factors, sampling at least 40,000 posterior
## samples is recommended.## Warning in summarize_brms(pa_obj_log_pressure, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 101.90*** | 6.39 | [ 89.74, 115.49] | 1.000 | [0.94, 1.07] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 7210 | 14382 |
| is_B | 137.60*** | 7.65 | [123.21, 153.70] | 1.000 | [0.94, 1.07] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 10736 | 18368 |
| is_A:pressure_self_cw | 0.99 | 0.05 | [ 0.87, 1.09] | 0.609 | [0.94, 1.07] | 0.764 | 0.006 | Very Strong Evidence for Null | 1.000 | 35375 | 26044 |
| is_B:pressure_self_cw | 0.99 | 0.04 | [ 0.90, 1.08] | 0.616 | [0.94, 1.07] | 0.843 | 0.005 | Very Strong Evidence for Null | 1.000 | 54399 | 31354 |
| is_A:pressure_partner_cw | 1.01 | 0.04 | [ 0.94, 1.09] | 0.603 | [0.94, 1.07] | 0.895 | 0.004 | Very Strong Evidence for Null | 1.000 | 56778 | 29304 |
| is_B:pressure_partner_cw | 1.02 | 0.06 | [ 0.92, 1.17] | 0.667 | [0.94, 1.07] | 0.709 | 0.006 | Very Strong Evidence for Null | 1.000 | 37773 | 26375 |
| is_A:pressure_self_cb | 0.84 | 0.28 | [ 0.43, 1.67] | 0.697 | [0.94, 1.07] | 0.135 | 0.017 | Very Strong Evidence for Null | 1.000 | 17538 | 23698 |
| is_B:pressure_self_cb | 0.94 | 0.22 | [ 0.58, 1.51] | 0.607 | [0.94, 1.07] | 0.210 | 0.012 | Very Strong Evidence for Null | 1.000 | 18333 | 23028 |
| is_A:pressure_partner_cb | 1.25 | 0.33 | [ 0.73, 2.14] | 0.797 | [0.94, 1.07] | 0.136 | 0.022 | Very Strong Evidence for Null | 1.000 | 16840 | 23213 |
| is_B:pressure_partner_cb | 1.20 | 0.36 | [ 0.65, 2.21] | 0.722 | [0.94, 1.07] | 0.139 | 0.014 | Very Strong Evidence for Null | 1.000 | 19355 | 24488 |
| is_A:day | 0.95 | 0.04 | [ 0.87, 1.05] | 0.842 | [0.94, 1.07] | 0.642 | 0.006 | Very Strong Evidence for Null | 1.000 | 83189 | 29167 |
| is_B:day | 0.92 | 0.04 | [ 0.84, 1.01] | 0.961 | [0.94, 1.07] | 0.360 | 0.018 | Very Strong Evidence for Null | 1.000 | 90539 | 28746 |
| is_A:weartime_self_cw | 1.00*** | 0.00 | [ 1.00, 1.00] | 1.000 | [0.94, 1.07] | 1.000 | 13.519 | Strong Evidence | 1.000 | 76844 | 32076 |
| is_B:weartime_self_cw | 1.00 | 0.00 | [ 1.00, 1.00] | 0.973 | [0.94, 1.07] | 1.000 | 0.025 | Very Strong Evidence for Null | 1.000 | 89204 | 27836 |
| is_A:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.631 | [0.94, 1.07] | 1.000 | 0.012 | Very Strong Evidence for Null | 1.000 | 9815 | 16497 |
| is_B:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.725 | [0.94, 1.07] | 1.000 | 0.014 | Very Strong Evidence for Null | 1.000 | 12846 | 21776 |
| sd(is_A) | 0.34 | 0.04 | [0.27, 0.44] | NA | NA | NA | NA | NA | 1.000 | 10546 | 18623 |
| sd(is_B) | 0.29 | 0.04 | [0.23, 0.38] | NA | NA | NA | NA | NA | 1.000 | 12764 | 19054 |
| sd(is_A:pressure_self_cw) | 0.08 | 0.07 | [0.00, 0.25] | NA | NA | NA | NA | NA | 1.000 | 14433 | 21148 |
| sd(is_B:pressure_self_cw) | 0.05 | 0.04 | [0.00, 0.17] | NA | NA | NA | NA | NA | 1.000 | 21197 | 20973 |
| sd(is_A:pressure_partner_cw) | 0.04 | 0.03 | [0.00, 0.13] | NA | NA | NA | NA | NA | 1.000 | 27025 | 21802 |
| sd(is_B:pressure_partner_cw) | 0.07 | 0.06 | [0.00, 0.27] | NA | NA | NA | NA | NA | 1.000 | 19755 | 22073 |
| sigma | 0.56 | 0.01 | [0.54, 0.57] | NA | NA | NA | NA | NA | 1.000 | 76294 | 29458 |
Pushing
formula <- bf(
pa_obj ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
barriers_self_cw + barriers_self_cb +
is_A:day + is_B:day +
is_A:weartime_self_cw + is_B:weartime_self_cw +
is_A:weartime_self_cb + is_B:weartime_self_cb +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = lognormal()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_obj_log_pushing <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = lognormal(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_pushing", suffix))
)## Warning: Rows containing NAs were excluded from the model.if (check_models) {
check_brms(pa_obj_log_pushing, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_obj_log_pushing, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 2.00 [1.88, 2.13] 1.41 0.50
## is_B 1.02 [1.00, 1.21] 1.01 0.98
## barriers_self_cw 1.07 [1.03, 1.13] 1.03 0.94
## barriers_self_cb 1.07 [1.04, 1.14] 1.03 0.93
## is_A:pushing_self_cw 1.05 [1.02, 1.12] 1.02 0.96
## is_B:pushing_self_cw 1.05 [1.02, 1.12] 1.02 0.96
## is_A:pushing_partner_cw 1.10 [1.07, 1.17] 1.05 0.91
## is_B:pushing_partner_cw 1.96 [1.85, 2.09] 1.40 0.51
## is_A:pushing_self_cb 1.83 [1.73, 1.94] 1.35 0.55
## is_B:pushing_self_cb 1.96 [1.85, 2.09] 1.40 0.51
## is_A:pushing_partner_cb 1.89 [1.78, 2.01] 1.37 0.53
## Tolerance 95% CI
## [0.47, 0.53]
## [0.83, 1.00]
## [0.88, 0.97]
## [0.88, 0.96]
## [0.89, 0.98]
## [0.89, 0.98]
## [0.86, 0.94]
## [0.48, 0.54]
## [0.51, 0.58]
## [0.48, 0.54]
## [0.50, 0.56]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 88%
## lognormal 9%
## weibull 3%
##
## Predicted Distribution of Response
##
## Distribution Probability
## lognormal 72%
## tweedie 16%
## neg. binomial (zero-infl.) 9%
##
## Divergences:## 0 of 40000 iterations ended with a divergence.##
## Tree depth:## 0 of 40000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 13, observations = 1594, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000925345 0.005646173
## sample estimates:
## outlier frequency (expected: 0.00331869510664994 )
## 0.008155583# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log_pushing$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_obj <- summarize_brms(
pa_obj_log_pushing,
stats_to_report = stats_to_report,
rope_range = rope_range_log,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.## Sampling priors, please wait...## Warning in summarize_brms(pa_obj_log_pushing, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 106.40*** | 7.11 | [ 93.02, 121.57] | 1.000 | [0.94, 1.06] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 15901 | 24332 |
| is_B | 146.91*** | 9.89 | [128.51, 168.01] | 1.000 | [0.94, 1.06] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 19659 | 25323 |
| barriers_self_cw | 0.95** | 0.02 | [ 0.91, 0.99] | 0.996 | [0.94, 1.06] | 0.647 | 0.177 | Moderate Evidence for Null | 1.000 | 91051 | 28693 |
| barriers_self_cb | 0.82* | 0.08 | [ 0.68, 0.99] | 0.980 | [0.94, 1.06] | 0.070 | 0.115 | Moderate Evidence for Null | 1.000 | 20765 | 25131 |
| is_A:pushing_self_cw | 1.07 | 0.04 | [ 0.99, 1.15] | 0.966 | [0.94, 1.06] | 0.433 | 0.046 | Strong Evidence for Null | 1.000 | 42743 | 30873 |
| is_B:pushing_self_cw | 1.00 | 0.03 | [ 0.94, 1.06] | 0.506 | [0.94, 1.06] | 0.951 | 0.007 | Very Strong Evidence for Null | 1.000 | 42212 | 30965 |
| is_A:pushing_partner_cw | 0.98 | 0.02 | [ 0.94, 1.04] | 0.728 | [0.94, 1.06] | 0.964 | 0.007 | Very Strong Evidence for Null | 1.000 | 67404 | 32270 |
| is_B:pushing_partner_cw | 1.06 | 0.04 | [ 0.99, 1.14] | 0.957 | [0.94, 1.06] | 0.491 | 0.035 | Strong Evidence for Null | 1.000 | 42448 | 31807 |
| is_A:pushing_self_cb | 0.98 | 0.28 | [ 0.55, 1.73] | 0.530 | [0.94, 1.06] | 0.174 | 0.015 | Very Strong Evidence for Null | 1.000 | 16613 | 23372 |
| is_B:pushing_self_cb | 0.98 | 0.31 | [ 0.52, 1.88] | 0.525 | [0.94, 1.06] | 0.154 | 0.018 | Very Strong Evidence for Null | 1.000 | 17222 | 23611 |
| is_A:pushing_partner_cb | 1.16 | 0.36 | [ 0.63, 2.18] | 0.688 | [0.94, 1.06] | 0.140 | 0.018 | Very Strong Evidence for Null | 1.000 | 15659 | 23690 |
| is_B:pushing_partner_cb | 1.42 | 0.41 | [ 0.79, 2.58] | 0.884 | [0.94, 1.06] | 0.080 | 0.033 | Strong Evidence for Null | 1.000 | 17956 | 24013 |
| is_A:day | 1.00 | 0.07 | [ 0.88, 1.14] | 0.523 | [0.94, 1.06] | 0.661 | 0.005 | Very Strong Evidence for Null | 1.000 | 81551 | 28339 |
| is_B:day | 0.89 | 0.06 | [ 0.78, 1.01] | 0.969 | [0.94, 1.06] | 0.186 | 0.031 | Strong Evidence for Null | 1.000 | 80514 | 29652 |
| is_A:weartime_self_cw | 1.00** | 0.00 | [ 1.00, 1.00] | 0.997 | [0.94, 1.06] | 1.000 | 0.237 | Moderate Evidence for Null | 1.000 | 74748 | 33237 |
| is_B:weartime_self_cw | 1.00* | 0.00 | [ 1.00, 1.00] | 0.993 | [0.94, 1.06] | 1.000 | 0.112 | Moderate Evidence for Null | 1.000 | 83919 | 30123 |
| is_A:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.703 | [0.94, 1.06] | 1.000 | 0.014 | Very Strong Evidence for Null | 1.000 | 16870 | 23704 |
| is_B:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.882 | [0.94, 1.06] | 1.000 | 0.028 | Very Strong Evidence for Null | 1.000 | 19627 | 25505 |
| sd(is_A) | 0.33 | 0.05 | [0.25, 0.44] | NA | NA | NA | NA | NA | 1.000 | 18055 | 25973 |
| sd(is_B) | 0.34 | 0.05 | [0.26, 0.45] | NA | NA | NA | NA | NA | 1.000 | 17389 | 24079 |
| sd(is_A:pushing_self_cw) | 0.10 | 0.05 | [0.01, 0.21] | NA | NA | NA | NA | NA | 1.000 | 9624 | 11024 |
| sd(is_B:pushing_self_cw) | 0.09 | 0.04 | [0.01, 0.17] | NA | NA | NA | NA | NA | 1.000 | 11822 | 9435 |
| sd(is_A:pushing_partner_cw) | 0.03 | 0.03 | [0.00, 0.10] | NA | NA | NA | NA | NA | 1.000 | 23687 | 23076 |
| sd(is_B:pushing_partner_cw) | 0.10 | 0.04 | [0.01, 0.18] | NA | NA | NA | NA | NA | 1.001 | 16684 | 13022 |
| sigma | 0.52 | 0.01 | [0.50, 0.54] | NA | NA | NA | NA | NA | 1.000 | 59441 | 29518 |
Affect
## [1] 0 5
Persuasion
formula <- bf(
aff ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b")
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = gaussian()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
mood_gauss_persuasion <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("mood_gauss_persuasion", suffix))
)## Warning: Rows containing NAs were excluded from the model.if (check_models) {
check_brms(mood_gauss_persuasion, log_pp_check = FALSE)
DHARMa.check_brms.all(mood_gauss_persuasion, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.17 [1.13, 1.21] 1.08 0.86
## is_B 1.04 [1.02, 1.09] 1.02 0.96
## is_A:persuasion_self_cw 1.06 [1.04, 1.11] 1.03 0.94
## is_B:persuasion_self_cw 1.07 [1.04, 1.11] 1.03 0.93
## is_A:persuasion_partner_cw 1.04 [1.02, 1.09] 1.02 0.96
## is_B:persuasion_partner_cw 2.85 [2.71, 2.99] 1.69 0.35
## is_A:persuasion_self_cb 2.80 [2.67, 2.94] 1.67 0.36
## Tolerance 95% CI
## [0.82, 0.88]
## [0.92, 0.98]
## [0.90, 0.96]
## [0.90, 0.96]
## [0.92, 0.98]
## [0.33, 0.37]
## [0.34, 0.38]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 47%
## normal 34%
## exponential 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## tweedie 41%
## beta-binomial 38%
## half-cauchy 6%
##
## Divergences:## 0 of 40000 iterations ended with a divergence.##
## Tree depth:## 0 of 40000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
##
## DHARMa outlier test based on exact binomial test with
## approximate expectations
##
## data: model.check
## outliers at both margin(s) = 39, observations = 3688, p-value =
## 2.231e-16
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.007530249 0.014428062
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.01057484summary_mood <- summarize_brms(
mood_gauss_persuasion,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = F) ## Model has no intercept. VIFs may not be sensible.## Sampling priors, please wait...| Est. | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 3.60*** | 0.12 | [ 3.35, 3.84] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.002 | 3494 | 6838 |
| is_B | 3.88*** | 0.14 | [ 3.60, 4.15] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 4288 | 9553 |
| is_A:persuasion_self_cw | 0.01 | 0.02 | [-0.04, 0.05] | 0.633 | [-0.11, 0.11] | 1.000 | 0.003 | Very Strong Evidence for Null | 1.000 | 53463 | 30917 |
| is_B:persuasion_self_cw | -0.01 | 0.02 | [-0.06, 0.03] | 0.717 | [-0.11, 0.11] | 1.000 | 0.003 | Very Strong Evidence for Null | 1.000 | 47802 | 29543 |
| is_A:persuasion_partner_cw | 0.05* | 0.02 | [ 0.00, 0.09] | 0.982 | [-0.11, 0.11] | 0.998 | 0.027 | Very Strong Evidence for Null | 1.000 | 45958 | 31085 |
| is_B:persuasion_partner_cw | 0.01 | 0.02 | [-0.04, 0.05] | 0.626 | [-0.11, 0.11] | 1.000 | 0.003 | Very Strong Evidence for Null | 1.000 | 55921 | 29328 |
| is_A:persuasion_self_cb | -0.29 | 0.33 | [-0.97, 0.38] | 0.810 | [-0.11, 0.11] | 0.185 | 0.023 | Very Strong Evidence for Null | 1.000 | 5518 | 10646 |
| is_B:persuasion_self_cb | 0.37 | 0.33 | [-0.30, 1.05] | 0.868 | [-0.11, 0.11] | 0.148 | 0.033 | Strong Evidence for Null | 1.001 | 6266 | 11337 |
| is_A:persuasion_partner_cb | 0.41 | 0.29 | [-0.18, 0.99] | 0.916 | [-0.11, 0.11] | 0.120 | 0.041 | Strong Evidence for Null | 1.000 | 5913 | 10851 |
| is_B:persuasion_partner_cb | -0.20 | 0.39 | [-0.98, 0.57] | 0.699 | [-0.11, 0.11] | 0.205 | 0.020 | Very Strong Evidence for Null | 1.000 | 6005 | 11939 |
| is_A:day | 0.45*** | 0.07 | [ 0.30, 0.58] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 51520 | 30176 |
| is_B:day | 0.08 | 0.07 | [-0.06, 0.22] | 0.871 | [-0.11, 0.11] | 0.673 | 0.006 | Very Strong Evidence for Null | 1.000 | 54047 | 28585 |
| sd(is_A) | 0.70 | 0.09 | [0.55, 0.90] | NA | NA | NA | NA | NA | 1.000 | 7045 | 14067 |
| sd(is_B) | 0.79 | 0.10 | [0.63, 1.03] | NA | NA | NA | NA | NA | 1.000 | 7522 | 12132 |
| sd(is_A:pushing_self_cw) | 0.11 | 0.06 | [0.01, 0.23] | NA | NA | NA | NA | NA | 1.000 | 10529 | 10208 |
| sd(is_B:pushing_self_cw) | 0.06 | 0.05 | [0.00, 0.18] | NA | NA | NA | NA | NA | 1.000 | 13285 | 16387 |
| sd(is_A:pushing_partner_cw) | 0.13 | 0.05 | [0.04, 0.24] | NA | NA | NA | NA | NA | 1.000 | 17381 | 12874 |
| sd(is_B:pushing_partner_cw) | 0.12 | 0.06 | [0.01, 0.25] | NA | NA | NA | NA | NA | 1.000 | 10329 | 7714 |
| sigma | 0.87 | 0.01 | [0.85, 0.89] | NA | NA | NA | NA | NA | 1.000 | 52078 | 28947 |
Pressure
formula <- bf(
aff ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b")
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = gaussian()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
mood_gauss_pressure <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("mood_gauss_pressure", suffix))
)## Warning: Rows containing NAs were excluded from the model.if (check_models) {
check_brms(mood_gauss_pressure, log_pp_check = FALSE)
DHARMa.check_brms.all(mood_gauss_pressure, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.19 [1.15, 1.23] 1.09 0.84
## is_B 1.01 [1.00, 1.21] 1.01 0.99
## is_A:pressure_self_cw 1.01 [1.00, 1.12] 1.01 0.99
## is_B:pressure_self_cw 1.01 [1.00, 1.12] 1.01 0.99
## is_A:pressure_partner_cw 1.01 [1.00, 1.42] 1.00 0.99
## Tolerance 95% CI
## [0.81, 0.87]
## [0.83, 1.00]
## [0.89, 1.00]
## [0.89, 1.00]
## [0.70, 1.00]
##
## Moderate Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_B:pressure_partner_cw 8.40 [7.94, 8.89] 2.90 0.12
## is_A:pressure_self_cb 8.05 [7.61, 8.52] 2.84 0.12
## Tolerance 95% CI
## [0.11, 0.13]
## [0.12, 0.13]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 47%
## normal 34%
## exponential 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## tweedie 41%
## beta-binomial 38%
## half-cauchy 6%
##
## Divergences:## 0 of 40000 iterations ended with a divergence.##
## Tree depth:## 0 of 40000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
##
## DHARMa outlier test based on exact binomial test with
## approximate expectations
##
## data: model.check
## outliers at both margin(s) = 41, observations = 3684, p-value <
## 2.2e-16
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.007998021 0.015068026
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.01112921summary_mood <- summarize_brms(
mood_gauss_pressure,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = F) ## Model has no intercept. VIFs may not be sensible.## Sampling priors, please wait...| Est. | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 3.64*** | 0.13 | [ 3.37, 3.89] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 2966 | 6814 |
| is_B | 3.88*** | 0.14 | [ 3.60, 4.17] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.003 | 3288 | 6776 |
| is_A:pressure_self_cw | -0.06 | 0.06 | [-0.17, 0.06] | 0.835 | [-0.11, 0.11] | 0.845 | 0.005 | Very Strong Evidence for Null | 1.000 | 58461 | 29843 |
| is_B:pressure_self_cw | 0.01 | 0.06 | [-0.10, 0.12] | 0.542 | [-0.11, 0.11] | 0.957 | 0.003 | Very Strong Evidence for Null | 1.000 | 44197 | 29836 |
| is_A:pressure_partner_cw | 0.05 | 0.06 | [-0.06, 0.16] | 0.808 | [-0.11, 0.11] | 0.864 | 0.005 | Very Strong Evidence for Null | 1.000 | 42613 | 30902 |
| is_B:pressure_partner_cw | 0.04 | 0.06 | [-0.07, 0.15] | 0.754 | [-0.11, 0.11] | 0.907 | 0.004 | Very Strong Evidence for Null | 1.000 | 68143 | 31795 |
| is_A:pressure_self_cb | -0.16 | 0.64 | [-1.48, 1.11] | 0.600 | [-0.11, 0.11] | 0.142 | 0.015 | Very Strong Evidence for Null | 1.000 | 10697 | 15995 |
| is_B:pressure_self_cb | 0.12 | 0.57 | [-1.03, 1.31] | 0.588 | [-0.11, 0.11] | 0.157 | 0.017 | Very Strong Evidence for Null | 1.000 | 10971 | 16363 |
| is_A:pressure_partner_cb | 0.22 | 0.50 | [-0.76, 1.27] | 0.674 | [-0.11, 0.11] | 0.164 | 0.016 | Very Strong Evidence for Null | 1.000 | 10831 | 15307 |
| is_B:pressure_partner_cb | -0.12 | 0.72 | [-1.61, 1.35] | 0.566 | [-0.11, 0.11] | 0.126 | 0.016 | Very Strong Evidence for Null | 1.000 | 11345 | 16597 |
| is_A:day | 0.41*** | 0.07 | [ 0.28, 0.55] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 49901 | 30936 |
| is_B:day | 0.09 | 0.07 | [-0.05, 0.22] | 0.897 | [-0.11, 0.11] | 0.647 | 0.006 | Very Strong Evidence for Null | 1.000 | 50514 | 29576 |
| sd(is_A) | 0.71 | 0.09 | [0.57, 0.93] | NA | NA | NA | NA | NA | 1.001 | 5363 | 10000 |
| sd(is_B) | 0.81 | 0.10 | [0.65, 1.06] | NA | NA | NA | NA | NA | 1.001 | 6354 | 13168 |
| sd(is_A:pushing_self_cw) | 0.12 | 0.06 | [0.01, 0.25] | NA | NA | NA | NA | NA | 1.001 | 9570 | 7098 |
| sd(is_B:pushing_self_cw) | 0.07 | 0.05 | [0.00, 0.18] | NA | NA | NA | NA | NA | 1.000 | 12893 | 14758 |
| sd(is_A:pushing_partner_cw) | 0.13 | 0.05 | [0.04, 0.24] | NA | NA | NA | NA | NA | 1.000 | 13798 | 9709 |
| sd(is_B:pushing_partner_cw) | 0.12 | 0.06 | [0.02, 0.25] | NA | NA | NA | NA | NA | 1.000 | 12399 | 9266 |
| sigma | 0.87 | 0.01 | [0.85, 0.89] | NA | NA | NA | NA | NA | 1.000 | 43311 | 29442 |
pushing
formula <- bf(
aff ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
barriers_self_cw + barriers_self_cb +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b")
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = gaussian()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
mood_gauss_pushing <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("mood_gauss_pushing", suffix))
)## Warning: Rows containing NAs were excluded from the model.if (check_models) {
check_brms(mood_gauss_pushing, log_pp_check = FALSE)
DHARMa.check_brms.all(mood_gauss_pushing, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.35 [1.29, 1.43] 1.16 0.74
## is_B 1.01 [1.00, 1.48] 1.01 0.99
## barriers_self_cw 1.09 [1.05, 1.15] 1.04 0.92
## barriers_self_cb 1.09 [1.05, 1.15] 1.04 0.92
## is_A:pushing_self_cw 1.05 [1.02, 1.12] 1.02 0.95
## is_B:pushing_self_cw 1.08 [1.05, 1.14] 1.04 0.93
## is_A:pushing_partner_cw 1.09 [1.05, 1.15] 1.04 0.92
## is_B:pushing_partner_cw 1.50 [1.42, 1.58] 1.22 0.67
## is_A:pushing_self_cb 1.54 [1.46, 1.62] 1.24 0.65
## Tolerance 95% CI
## [0.70, 0.77]
## [0.68, 1.00]
## [0.87, 0.95]
## [0.87, 0.95]
## [0.89, 0.98]
## [0.88, 0.96]
## [0.87, 0.95]
## [0.63, 0.70]
## [0.62, 0.68]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 41%
## normal 41%
## exponential 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## tweedie 41%
## beta-binomial 38%
## half-cauchy 6%
##
## Divergences:## 0 of 40000 iterations ended with a divergence.##
## Tree depth:## 0 of 40000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
##
## DHARMa outlier test based on exact binomial test with
## approximate expectations
##
## data: model.check
## outliers at both margin(s) = 11, observations = 1776, p-value =
## 0.00112
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.003095804 0.011055146
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.006193694summary_mood <- summarize_brms(
mood_gauss_pushing,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = F) ## Model has no intercept. VIFs may not be sensible.## Sampling priors, please wait...| Est. | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 3.66*** | 0.12 | [ 3.41, 3.90] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 7565 | 14933 |
| is_B | 3.84*** | 0.14 | [ 3.55, 4.13] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 7731 | 13693 |
| barriers_self_cw | -0.29*** | 0.03 | [-0.35, -0.23] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 78719 | 29655 |
| barriers_self_cb | -0.50* | 0.22 | [-0.95, -0.05] | 0.984 | [-0.11, 0.11] | 0.039 | 0.167 | Moderate Evidence for Null | 1.000 | 8107 | 14716 |
| is_A:pushing_self_cw | 0.09 | 0.05 | [ 0.00, 0.18] | 0.969 | [-0.11, 0.11] | 0.659 | 0.034 | Strong Evidence for Null | 1.000 | 42568 | 31360 |
| is_B:pushing_self_cw | -0.01 | 0.04 | [-0.09, 0.06] | 0.640 | [-0.11, 0.11] | 0.990 | 0.005 | Very Strong Evidence for Null | 1.000 | 48493 | 32295 |
| is_A:pushing_partner_cw | 0.09 | 0.04 | [ 0.00, 0.18] | 0.974 | [-0.11, 0.11] | 0.678 | 0.039 | Strong Evidence for Null | 1.000 | 32188 | 29785 |
| is_B:pushing_partner_cw | 0.13* | 0.05 | [ 0.02, 0.23] | 0.985 | [-0.11, 0.11] | 0.376 | 0.077 | Strong Evidence for Null | 1.000 | 33217 | 29086 |
| is_A:pushing_self_cb | -0.23 | 0.48 | [-1.20, 0.73] | 0.681 | [-0.11, 0.11] | 0.162 | 0.015 | Very Strong Evidence for Null | 1.000 | 11236 | 17666 |
| is_B:pushing_self_cb | 0.90 | 0.67 | [-0.47, 2.27] | 0.905 | [-0.11, 0.11] | 0.053 | 0.047 | Strong Evidence for Null | 1.000 | 9401 | 14271 |
| is_A:pushing_partner_cb | 0.77 | 0.52 | [-0.29, 1.83] | 0.927 | [-0.11, 0.11] | 0.056 | 0.042 | Strong Evidence for Null | 1.001 | 9141 | 15831 |
| is_B:pushing_partner_cb | -0.32 | 0.63 | [-1.57, 0.96] | 0.693 | [-0.11, 0.11] | 0.124 | 0.022 | Very Strong Evidence for Null | 1.000 | 11283 | 17409 |
| is_A:day | 0.34*** | 0.10 | [ 0.15, 0.54] | 1.000 | [-0.11, 0.11] | 0.008 | 1.941 | Weak Evidence | 1.000 | 70757 | 30332 |
| is_B:day | 0.08 | 0.10 | [-0.11, 0.27] | 0.793 | [-0.11, 0.11] | 0.609 | 0.006 | Very Strong Evidence for Null | 1.000 | 69006 | 28830 |
| sd(is_A) | 0.64 | 0.09 | [0.50, 0.84] | NA | NA | NA | NA | NA | 1.001 | 10062 | 18475 |
| sd(is_B) | 0.79 | 0.10 | [0.63, 1.04] | NA | NA | NA | NA | NA | 1.000 | 10080 | 17533 |
| sd(is_A:pushing_self_cw) | 0.09 | 0.07 | [0.01, 0.24] | NA | NA | NA | NA | NA | 1.000 | 12455 | 15086 |
| sd(is_B:pushing_self_cw) | 0.07 | 0.05 | [0.00, 0.20] | NA | NA | NA | NA | NA | 1.000 | 14139 | 16967 |
| sd(is_A:pushing_partner_cw) | 0.12 | 0.06 | [0.01, 0.24] | NA | NA | NA | NA | NA | 1.000 | 13704 | 10047 |
| sd(is_B:pushing_partner_cw) | 0.14 | 0.07 | [0.02, 0.29] | NA | NA | NA | NA | NA | 1.000 | 14270 | 12569 |
| sigma | 0.83 | 0.01 | [0.80, 0.86] | NA | NA | NA | NA | NA | 1.000 | 60197 | 29169 |
Reactance
## [1] 0 5
introduce_binary_reactance <- function(data) {
data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
return(data)
}
df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
for (i in seq_along(implist)) {
implist[[i]] <- introduce_binary_reactance(implist[[i]])
}
}Persuasion
formula <- bf(
is_reactance ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = bernoulli()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
is_reactance_persuasion <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::bernoulli(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("is_reactance_persuasion", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.if (check_models) {
check_brms(is_reactance_persuasion)
DHARMa.check_brms.all(is_reactance_persuasion, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.82 [1.71, 1.96] 1.35 0.55
## is_B 1.39 [1.32, 1.49] 1.18 0.72
## is_A:persuasion_self_cw 1.15 [1.10, 1.23] 1.07 0.87
## is_B:persuasion_self_cw 1.04 [1.01, 1.14] 1.02 0.96
## is_A:persuasion_partner_cw 1.06 [1.02, 1.14] 1.03 0.95
## is_B:persuasion_partner_cw 1.51 [1.43, 1.62] 1.23 0.66
## is_A:persuasion_self_cb 1.78 [1.67, 1.91] 1.33 0.56
## Tolerance 95% CI
## [0.51, 0.58]
## [0.67, 0.76]
## [0.82, 0.91]
## [0.88, 0.99]
## [0.88, 0.98]
## [0.62, 0.70]
## [0.52, 0.60]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## normal 44%
## cauchy 12%
## gamma 9%
##
## Predicted Distribution of Response
##
## Distribution Probability
## bernoulli 81%
## beta-binomial 12%
## binomial 6%
##
## Divergences:## 0 of 40000 iterations ended with a divergence.##
## Tree depth:## 0 of 40000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 5 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
##
## DHARMa outlier test based on exact binomial test with
## approximate expectations
##
## data: model.check
## outliers at both margin(s) = 1, observations = 756, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.0000334886 0.0073476538
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.001322751summary_is_reactance <- summarize_brms(
is_reactance_persuasion,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.## Sampling priors, please wait...| OR | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.18*** | 0.08 | [0.07, 0.43] | 1.000 | [0.83, 1.20] | 0.000 | 36.028 | Very Strong Evidence | 1.000 | 24897 | 29836 |
| is_B | 0.30** | 0.12 | [0.12, 0.65] | 0.999 | [0.83, 1.20] | 0.004 | 2.802 | Weak Evidence | 1.000 | 29961 | 29498 |
| is_A:persuasion_self_cw | 0.90 | 0.13 | [0.64, 1.18] | 0.773 | [0.83, 1.20] | 0.669 | 0.071 | Strong Evidence for Null | 1.000 | 22239 | 26651 |
| is_B:persuasion_self_cw | 0.78 | 0.12 | [0.54, 1.02] | 0.965 | [0.83, 1.20] | 0.310 | 0.293 | Moderate Evidence for Null | 1.000 | 23382 | 19058 |
| is_A:persuasion_partner_cw | 0.89 | 0.16 | [0.58, 1.29] | 0.740 | [0.83, 1.20] | 0.579 | 0.077 | Strong Evidence for Null | 1.000 | 30521 | 22671 |
| is_B:persuasion_partner_cw | 1.27 | 0.24 | [0.86, 1.97] | 0.897 | [0.83, 1.20] | 0.356 | 0.165 | Moderate Evidence for Null | 1.000 | 25245 | 25126 |
| is_A:persuasion_self_cb | 5.13 | 4.34 | [0.96, 35.39] | 0.972 | [0.83, 1.20] | 0.025 | 0.569 | Weak Evidence for Null | 1.000 | 18611 | 24445 |
| is_B:persuasion_self_cb | 2.76 | 2.14 | [0.62, 15.25] | 0.910 | [0.83, 1.20] | 0.078 | 0.216 | Moderate Evidence for Null | 1.000 | 20515 | 24894 |
| is_A:persuasion_partner_cb | 1.58 | 1.13 | [0.40, 7.92] | 0.746 | [0.83, 1.20] | 0.169 | 0.114 | Moderate Evidence for Null | 1.000 | 16564 | 19332 |
| is_B:persuasion_partner_cb | 2.13 | 1.94 | [0.34, 14.06] | 0.798 | [0.83, 1.20] | 0.109 | 0.122 | Moderate Evidence for Null | 1.000 | 18086 | 24106 |
| is_A:day | 2.57 | 1.36 | [0.91, 7.43] | 0.962 | [0.83, 1.20] | 0.058 | 0.194 | Moderate Evidence for Null | 1.000 | 56235 | 31190 |
| is_B:day | 1.04 | 0.59 | [0.33, 3.17] | 0.523 | [0.83, 1.20] | 0.248 | 0.045 | Strong Evidence for Null | 1.000 | 52331 | 30861 |
| sd(is_A) | 1.38 | 0.42 | [0.68, 2.40] | NA | NA | NA | NA | NA | 1.000 | 14059 | 19287 |
| sd(is_B) | 1.36 | 0.33 | [0.82, 2.18] | NA | NA | NA | NA | NA | 1.000 | 16662 | 23957 |
| sd(is_A:persuasion_self_cw) | 0.31 | 0.18 | [0.03, 0.73] | NA | NA | NA | NA | NA | 1.001 | 9510 | 11805 |
| sd(is_B:persuasion_self_cw) | 0.35 | 0.26 | [0.02, 0.95] | NA | NA | NA | NA | NA | 1.000 | 5934 | 13700 |
| sd(is_A:persuasion_partner_cw) | 0.47 | 0.27 | [0.05, 1.19] | NA | NA | NA | NA | NA | 1.000 | 12291 | 11597 |
| sd(is_B:persuasion_partner_cw) | 0.65 | 0.29 | [0.19, 1.40] | NA | NA | NA | NA | NA | 1.000 | 16040 | 15521 |
Pressure
formula <- bf(
is_reactance ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = bernoulli()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
is_reactance_pressure <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::bernoulli(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("is_reactance_pressure", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.if (check_models) {
check_brms(is_reactance_pressure)
DHARMa.check_brms.all(is_reactance_pressure, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.83 [1.72, 1.96] 1.35 0.55
## is_B 1.02 [1.00, 1.23] 1.01 0.98
## is_A:pressure_self_cw 1.01 [1.00, 1.59] 1.01 0.99
## is_B:pressure_self_cw 1.01 [1.00, 2.01] 1.01 0.99
## is_A:pressure_partner_cw 1.01 [1.00, 36.14] 1.00 0.99
## is_B:pressure_partner_cw 2.51 [2.33, 2.71] 1.58 0.40
## is_A:pressure_self_cb 2.97 [2.76, 3.22] 1.72 0.34
## Tolerance 95% CI
## [0.51, 0.58]
## [0.82, 1.00]
## [0.63, 1.00]
## [0.50, 1.00]
## [0.03, 1.00]
## [0.37, 0.43]
## [0.31, 0.36]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## normal 34%
## beta 16%
## cauchy 12%
##
## Predicted Distribution of Response
##
## Distribution Probability
## bernoulli 81%
## beta-binomial 12%
## binomial 6%
##
## Divergences:## 0 of 40000 iterations ended with a divergence.##
## Tree depth:## 0 of 40000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 23 observations with a pareto_k > 0.7 in model 'model'.
## We recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
##
## DHARMa outlier test based on exact binomial test with
## approximate expectations
##
## data: model.check
## outliers at both margin(s) = 1, observations = 756, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.0000334886 0.0073476538
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.001322751summary_is_reactance <- summarize_brms(
is_reactance_pressure,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.## Sampling priors, please wait...| OR | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.15*** | 0.06 | [0.06, 0.32] | 1.000 | [0.83, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 27201 | 29563 |
| is_B | 0.21*** | 0.07 | [0.10, 0.40] | 1.000 | [0.83, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 33935 | 30484 |
| is_A:pressure_self_cw | 2.48 | 1.27 | [0.82, 11.13] | 0.957 | [0.83, 1.20] | 0.046 | 0.724 | Weak Evidence for Null | 1.000 | 21167 | 16301 |
| is_B:pressure_self_cw | 1.84 | 0.67 | [0.78, 4.37] | 0.935 | [0.83, 1.20] | 0.096 | 0.359 | Weak Evidence for Null | 1.000 | 27629 | 24426 |
| is_A:pressure_partner_cw | 0.87 | 0.65 | [0.09, 3.62] | 0.581 | [0.83, 1.20] | 0.205 | 0.128 | Moderate Evidence for Null | 1.000 | 24646 | 22567 |
| is_B:pressure_partner_cw | 1.56 | 0.77 | [0.39, 5.39] | 0.805 | [0.83, 1.20] | 0.170 | 0.128 | Moderate Evidence for Null | 1.000 | 30905 | 28417 |
| is_A:pressure_self_cb | 16.91 | 28.48 | [0.78, 1025.11] | 0.965 | [0.83, 1.20] | 0.018 | 0.537 | Weak Evidence for Null | 1.000 | 26287 | 20773 |
| is_B:pressure_self_cb | 23.56* | 34.63 | [1.93, 874.04] | 0.994 | [0.83, 1.20] | 0.006 | 1.817 | Weak Evidence | 1.000 | 25130 | 23908 |
| is_A:pressure_partner_cb | 1.05 | 1.44 | [0.05, 17.51] | 0.513 | [0.83, 1.20] | 0.105 | 0.120 | Moderate Evidence for Null | 1.000 | 24236 | 22816 |
| is_B:pressure_partner_cb | 0.36 | 0.65 | [0.00, 9.70] | 0.725 | [0.83, 1.20] | 0.071 | 0.105 | Moderate Evidence for Null | 1.000 | 29760 | 25453 |
| is_A:day | 2.24 | 1.18 | [0.79, 6.37] | 0.936 | [0.83, 1.20] | 0.087 | 0.138 | Moderate Evidence for Null | 1.000 | 59228 | 29462 |
| is_B:day | 1.34 | 0.70 | [0.48, 3.72] | 0.710 | [0.83, 1.20] | 0.234 | 0.050 | Strong Evidence for Null | 1.000 | 64117 | 29915 |
| sd(is_A) | 1.46 | 0.39 | [0.84, 2.44] | NA | NA | NA | NA | NA | 1.000 | 14625 | 22279 |
| sd(is_B) | 1.14 | 0.28 | [0.67, 1.83] | NA | NA | NA | NA | NA | 1.000 | 15965 | 24630 |
| sd(is_A:pressure_self_cw) | 1.30 | 1.07 | [0.06, 3.91] | NA | NA | NA | NA | NA | 1.000 | 9827 | 16955 |
| sd(is_B:pressure_self_cw) | 1.03 | 0.55 | [0.17, 2.56] | NA | NA | NA | NA | NA | 1.000 | 13715 | 12493 |
| sd(is_A:pressure_partner_cw) | 1.81 | 1.16 | [0.13, 4.53] | NA | NA | NA | NA | NA | 1.000 | 13418 | 11212 |
| sd(is_B:pressure_partner_cw) | 0.74 | 0.66 | [0.03, 2.89] | NA | NA | NA | NA | NA | 1.000 | 22546 | 21566 |
pushing
formula <- bf(
is_reactance ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
barriers_self_cw + barriers_self_cb +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = bernoulli()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
is_reactance_pushing <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::bernoulli(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("is_reactance_pushing", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.if (check_models) {
check_brms(is_reactance_pushing)
DHARMa.check_brms.all(is_reactance_pushing, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.85 [1.75, 1.97] 1.36 0.54
## is_B 1.06 [1.03, 1.13] 1.03 0.94
## barriers_self_cw 1.08 [1.05, 1.15] 1.04 0.92
## barriers_self_cb 1.24 [1.19, 1.31] 1.11 0.81
## is_A:pushing_self_cw 1.14 [1.10, 1.21] 1.07 0.87
## is_B:pushing_self_cw 1.01 [1.00, 1.89] 1.00 0.99
## is_A:pushing_partner_cw 1.05 [1.02, 1.12] 1.02 0.96
## is_B:pushing_partner_cw 1.53 [1.45, 1.62] 1.24 0.65
## is_A:pushing_self_cb 1.33 [1.27, 1.40] 1.15 0.75
## Tolerance 95% CI
## [0.51, 0.57]
## [0.89, 0.97]
## [0.87, 0.96]
## [0.77, 0.84]
## [0.83, 0.91]
## [0.53, 1.00]
## [0.89, 0.98]
## [0.62, 0.69]
## [0.71, 0.79]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## normal 44%
## cauchy 12%
## tweedie 9%
##
## Predicted Distribution of Response
##
## Distribution Probability
## bernoulli 75%
## beta-binomial 16%
## binomial 9%
##
## Divergences:## 0 of 40000 iterations ended with a divergence.##
## Tree depth:## 0 of 40000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 3 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
##
## DHARMa outlier test based on exact binomial test with
## approximate expectations
##
## data: model.check
## outliers at both margin(s) = 0, observations = 486, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.000000000 0.007561553
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0summary_is_reactance <- summarize_brms(
is_reactance_pushing,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.## Sampling priors, please wait...| OR | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.13*** | 0.07 | [0.04, 0.33] | 1.000 | [0.83, 1.20] | 0.000 | 75.098 | Very Strong Evidence | 1.000 | 21167 | 22208 |
| is_B | 0.16*** | 0.07 | [0.06, 0.35] | 1.000 | [0.83, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 28688 | 28514 |
| barriers_self_cw | 0.80 | 0.17 | [0.53, 1.19] | 0.865 | [0.83, 1.20] | 0.388 | 0.098 | Strong Evidence for Null | 1.000 | 43333 | 29118 |
| barriers_self_cb | 2.74 | 2.06 | [0.61, 12.70] | 0.910 | [0.83, 1.20] | 0.077 | 0.221 | Moderate Evidence for Null | 1.000 | 20027 | 24196 |
| is_A:pushing_self_cw | 1.35 | 0.25 | [0.93, 2.00] | 0.944 | [0.83, 1.20] | 0.258 | 0.217 | Moderate Evidence for Null | 1.000 | 35407 | 26611 |
| is_B:pushing_self_cw | 1.44 | 0.29 | [0.96, 2.23] | 0.964 | [0.83, 1.20] | 0.167 | 0.440 | Weak Evidence for Null | 1.000 | 25829 | 25680 |
| is_A:pushing_partner_cw | 0.77 | 0.24 | [0.32, 1.36] | 0.815 | [0.83, 1.20] | 0.338 | 0.129 | Moderate Evidence for Null | 1.000 | 19974 | 15679 |
| is_B:pushing_partner_cw | 1.15 | 0.23 | [0.75, 1.74] | 0.752 | [0.83, 1.20] | 0.517 | 0.079 | Strong Evidence for Null | 1.000 | 35235 | 26529 |
| is_A:pushing_self_cb | 38.45* | 65.62 | [2.09, 2544.92] | 0.993 | [0.83, 1.20] | 0.005 | 2.260 | Weak Evidence | 1.000 | 12493 | 16241 |
| is_B:pushing_self_cb | 2.37 | 3.65 | [0.13, 75.77] | 0.717 | [0.83, 1.20] | 0.082 | 0.123 | Moderate Evidence for Null | 1.000 | 17952 | 21560 |
| is_A:pushing_partner_cb | 0.08 | 0.16 | [0.00, 4.46] | 0.898 | [0.83, 1.20] | 0.032 | 0.268 | Moderate Evidence for Null | 1.000 | 16536 | 21879 |
| is_B:pushing_partner_cb | 1.25 | 1.87 | [0.05, 26.76] | 0.561 | [0.83, 1.20] | 0.097 | 0.100 | Strong Evidence for Null | 1.000 | 19752 | 22818 |
| is_A:day | 1.28 | 0.94 | [0.29, 5.36] | 0.633 | [0.83, 1.20] | 0.184 | 0.059 | Strong Evidence for Null | 1.000 | 39476 | 29465 |
| is_B:day | 1.40 | 0.92 | [0.39, 5.18] | 0.696 | [0.83, 1.20] | 0.194 | 0.064 | Strong Evidence for Null | 1.000 | 47667 | 31034 |
| sd(is_A) | 1.52 | 0.52 | [0.61, 2.78] | NA | NA | NA | NA | NA | 1.000 | 8591 | 8207 |
| sd(is_B) | 1.26 | 0.38 | [0.61, 2.16] | NA | NA | NA | NA | NA | 1.000 | 14530 | 19907 |
| sd(is_A:pushing_self_cw) | 0.26 | 0.23 | [0.01, 0.91] | NA | NA | NA | NA | NA | 1.000 | 10496 | 17794 |
| sd(is_B:pushing_self_cw) | 0.56 | 0.30 | [0.07, 1.32] | NA | NA | NA | NA | NA | 1.001 | 8769 | 8991 |
| sd(is_A:pushing_partner_cw) | 0.64 | 0.46 | [0.04, 1.94] | NA | NA | NA | NA | NA | 1.000 | 10726 | 13385 |
| sd(is_B:pushing_partner_cw) | 0.28 | 0.25 | [0.01, 1.05] | NA | NA | NA | NA | NA | 1.000 | 14535 | 18263 |
Analyses were conducted using the R Statistical language (version 4.4.2; R Core Team, 2024) on Windows 11 x64 (build 26100)
- rmarkdown (version 2.29; Allaire J et al., 2024)
- beepr (version 2.0; Bååth R, 2024)
- R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
- R.oo (version 1.27.0; Bengtsson H, 2003)
- R.utils (version 2.12.3; Bengtsson H, 2023)
- brms (version 2.22.0; Bürkner P, 2017)
- posterior (version 1.6.0; Bürkner P et al., 2024)
- extrafont (version 0.19; Chang W, 2023)
- digest (version 0.6.37; Eddelbuettel D, 2024)
- Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
- bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
- lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
- DHARMa (version 0.4.7; Hartig F, 2024)
- rlang (version 1.1.4; Henry L, Wickham H, 2024)
- tidybayes (version 3.0.7; Kay M, 2024)
- wbCorr (version 0.1.22; Küng P, 2023)
- see (version 0.9.0; Lüdecke D et al., 2021)
- tibble (version 3.2.1; Müller K, Wickham H, 2023)
- patchwork (version 1.3.0; Pedersen T, 2024)
- R (version 4.4.2; R Core Team, 2024)
- openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
- plotly (version 4.10.4; Sievert C, 2020)
- ggplot2 (version 3.5.1; Wickham H, 2016)
- forcats (version 1.0.0; Wickham H, 2023)
- stringr (version 1.5.1; Wickham H, 2023)
- rvest (version 1.0.4; Wickham H, 2024)
- tidyverse (version 2.0.0; Wickham H et al., 2019)
- readxl (version 1.4.3; Wickham H, Bryan J, 2023)
- dplyr (version 1.1.4; Wickham H et al., 2023)
- purrr (version 1.0.2; Wickham H, Henry L, 2023)
- readr (version 2.1.5; Wickham H et al., 2024)
- xml2 (version 1.3.6; Wickham H et al., 2023)
- tidyr (version 1.3.1; Wickham H et al., 2024)
- ggridges (version 0.5.6; Wilke C, 2024)
- knitr (version 1.49; Xie Y, 2024)
- kableExtra (version 1.4.0; Zhu H, 2024)
- Allaire J, Xie Y, Dervieux C, McPherson J, Luraschi J, Ushey K, Atkins A, Wickham H, Cheng J, Chang W, Iannone R (2024). rmarkdown: Dynamic Documents for R. R package version 2.29, https://github.com/rstudio/rmarkdown. Xie Y, Allaire J, Grolemund G (2018). R Markdown: The Definitive Guide. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9781138359338, https://bookdown.org/yihui/rmarkdown. Xie Y, Dervieux C, Riederer E (2020). R Markdown Cookbook. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9780367563837, https://bookdown.org/yihui/rmarkdown-cookbook.
- Bååth R (2024). beepr: Easily Play Notification Sounds on any Platform. R package version 2.0, https://github.com/rasmusab/beepr.
- Bengtsson H (2003). “The R.oo package - Object-Oriented Programming with References Using Standard R Code.” In Hornik K, Leisch F, Zeileis A (eds.), Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003). https://www.r-project.org/conferences/DSC-2003/Proceedings/Bengtsson.pdf.
- Bengtsson H (2003). “The R.oo package - Object-Oriented Programming with References Using Standard R Code.” In Hornik K, Leisch F, Zeileis A (eds.), Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003). https://www.r-project.org/conferences/DSC-2003/Proceedings/Bengtsson.pdf.
- Bengtsson H (2023). R.utils: Various Programming Utilities. R package version 2.12.3, https://github.com/HenrikBengtsson/R.utils, https://henrikbengtsson.github.io/R.utils/.
- Bürkner P (2017). “brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software, 80(1), 1-28. doi:10.18637/jss.v080.i01 https://doi.org/10.18637/jss.v080.i01. Bürkner P (2018). “Advanced Bayesian Multilevel Modeling with the R Package brms.” The R Journal, 10(1), 395-411. doi:10.32614/RJ-2018-017 https://doi.org/10.32614/RJ-2018-017. Bürkner P (2021). “Bayesian Item Response Modeling in R with brms and Stan.” Journal of Statistical Software, 100(5), 1-54. doi:10.18637/jss.v100.i05 https://doi.org/10.18637/jss.v100.i05.
- Bürkner P, Gabry J, Kay M, Vehtari A (2024). “posterior: Tools for Working with Posterior Distributions.” R package version 1.6.0, https://mc-stan.org/posterior/. Vehtari A, Gelman A, Simpson D, Carpenter B, Bürkner P (2021). “Rank-normalization, folding, and localization: An improved Rhat for assessing convergence of MCMC (with discussion).” Bayesian Analysis.
- Chang W (2023). extrafont: Tools for Using Fonts. R package version 0.19, https://github.com/wch/extrafont.
- Eddelbuettel D (2024). digest: Create Compact Hash Digests of R Objects. R package version 0.6.37, https://dirk.eddelbuettel.com/code/digest.html, https://github.com/eddelbuettel/digest.
- Eddelbuettel D, Francois R, Allaire J, Ushey K, Kou Q, Russell N, Ucar I, Bates D, Chambers J (2024). Rcpp: Seamless R and C++ Integration. R package version 1.0.13-1, https://dirk.eddelbuettel.com/code/rcpp.html, https://github.com/RcppCore/Rcpp, https://www.rcpp.org. Eddelbuettel D, François R (2011). “Rcpp: Seamless R and C++ Integration.” Journal of Statistical Software, 40(8), 1-18. doi:10.18637/jss.v040.i08 https://doi.org/10.18637/jss.v040.i08. Eddelbuettel D (2013). Seamless R and C++ Integration with Rcpp. Springer, New York. doi:10.1007/978-1-4614-6868-4 https://doi.org/10.1007/978-1-4614-6868-4, ISBN 978-1-4614-6867-7. Eddelbuettel D, Balamuta J (2018). “Extending R with C++: A Brief Introduction to Rcpp.” The American Statistician, 72(1), 28-36. doi:10.1080/00031305.2017.1375990 https://doi.org/10.1080/00031305.2017.1375990.
- Gabry J, Mahr T (2024). “bayesplot: Plotting for Bayesian Models.” R package version 1.11.1, https://mc-stan.org/bayesplot/. Gabry J, Simpson D, Vehtari A, Betancourt M, Gelman A (2019). “Visualization in Bayesian workflow.” J. R. Stat. Soc. A, 182, 389-402. doi:10.1111/rssa.12378 https://doi.org/10.1111/rssa.12378.
- Grolemund G, Wickham H (2011). “Dates and Times Made Easy with lubridate.” Journal of Statistical Software, 40(3), 1-25. https://www.jstatsoft.org/v40/i03/.
- Hartig F (2024). DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models. R package version 0.4.7, http://florianhartig.github.io/DHARMa/.
- Henry L, Wickham H (2024). rlang: Functions for Base Types and Core R and ‘Tidyverse’ Features. R package version 1.1.4, https://github.com/r-lib/rlang, https://rlang.r-lib.org.
- Kay M (2024). tidybayes: Tidy Data and Geoms for Bayesian Models. doi:10.5281/zenodo.1308151 https://doi.org/10.5281/zenodo.1308151, R package version 3.0.7, http://mjskay.github.io/tidybayes/.
- Küng P (2023). wbCorr: Bivariate Within- and Between-Cluster Correlations. University of Zürich. R package version 0.1.22. https://github.com/Pascal-Kueng/wbCorr.
- Lüdecke D, Patil I, Ben-Shachar M, Wiernik B, Waggoner P, Makowski D (2021). “see: An R Package for Visualizing Statistical Models.” Journal of Open Source Software, 6(64), 3393. doi:10.21105/joss.03393 https://doi.org/10.21105/joss.03393.
- Müller K, Wickham H (2023). tibble: Simple Data Frames. R package version 3.2.1, https://github.com/tidyverse/tibble, https://tibble.tidyverse.org/.
- Pedersen T (2024). patchwork: The Composer of Plots. R package version 1.3.0, https://github.com/thomasp85/patchwork, https://patchwork.data-imaginist.com.
- R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.
- Schauberger P, Walker A (2024). openxlsx: Read, Write and Edit xlsx Files. R package version 4.2.7.1, https://github.com/ycphs/openxlsx, https://ycphs.github.io/openxlsx/index.html.
- Sievert C (2020). Interactive Web-Based Data Visualization with R, plotly, and shiny. Chapman and Hall/CRC. ISBN 9781138331457, https://plotly-r.com.
- Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, https://ggplot2.tidyverse.org.
- Wickham H (2023). forcats: Tools for Working with Categorical Variables (Factors). R package version 1.0.0, https://github.com/tidyverse/forcats, https://forcats.tidyverse.org/.
- Wickham H (2023). stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.5.1, https://github.com/tidyverse/stringr, https://stringr.tidyverse.org.
- Wickham H (2024). rvest: Easily Harvest (Scrape) Web Pages. R package version 1.0.4, https://github.com/tidyverse/rvest, https://rvest.tidyverse.org/.
- Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). “Welcome to the tidyverse.” Journal of Open Source Software, 4(43), 1686. doi:10.21105/joss.01686 https://doi.org/10.21105/joss.01686.
- Wickham H, Bryan J (2023). readxl: Read Excel Files. R package version 1.4.3, https://github.com/tidyverse/readxl, https://readxl.tidyverse.org.
- Wickham H, François R, Henry L, Müller K, Vaughan D (2023). dplyr: A Grammar of Data Manipulation. R package version 1.1.4, https://github.com/tidyverse/dplyr, https://dplyr.tidyverse.org.
- Wickham H, Henry L (2023). purrr: Functional Programming Tools. R package version 1.0.2, https://github.com/tidyverse/purrr, https://purrr.tidyverse.org/.
- Wickham H, Hester J, Bryan J (2024). readr: Read Rectangular Text Data. R package version 2.1.5, https://github.com/tidyverse/readr, https://readr.tidyverse.org.
- Wickham H, Hester J, Ooms J (2023). xml2: Parse XML. R package version 1.3.6, https://github.com/r-lib/xml2, https://xml2.r-lib.org/.
- Wickham H, Vaughan D, Girlich M (2024). tidyr: Tidy Messy Data. R package version 1.3.1, https://github.com/tidyverse/tidyr, https://tidyr.tidyverse.org.
- Wilke C (2024). ggridges: Ridgeline Plots in ‘ggplot2’. R package version 0.5.6, https://wilkelab.org/ggridges/.
- Xie Y (2024). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.49, https://yihui.org/knitr/. Xie Y (2015). Dynamic Documents with R and knitr, 2nd edition. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, https://yihui.org/knitr/. Xie Y (2014). “knitr: A Comprehensive Tool for Reproducible Research in R.” In Stodden V, Leisch F, Peng RD (eds.), Implementing Reproducible Computational Research. Chapman and Hall/CRC. ISBN 978-1466561595.
- Zhu H (2024). kableExtra: Construct Complex Table with ‘kable’ and Pipe Syntax. R package version 1.4.0, https://github.com/haozhu233/kableExtra, http://haozhu233.github.io/kableExtra/.